# Figure 3 in manuscript

# Shows the xylem vulnerability curves for six species of subtropical thicket

# 29 June 2023
setwd("C://Users/rober/OneDrive/Documents/Thicket/raw data/Figures/")
jpeg(filename = "Figure3.jpg",
     width = 7, height = 10, units = "cm", pointsize = 1,
     bg = "white", res = 800, family = "", restoreConsole = TRUE,
     type = c("windows"),antialias = "cleartype")
# tiff(filename = "Figure2.tif",
#      width = 6, height = 10, units = "cm", pointsize = 1,
#      compression = c("lzw"),
#      bg = "white", res = 800, family = "", restoreConsole = TRUE,
#      type = c("windows"),antialias = "cleartype")
cols <- c("#F2EDD8","#F2EDD8","#F2EDD8","#F2EDD8","#F2EDD8","#F2EDD8")
par(mfrow=c(3,2),oma=c(6,5,0.15,0.15),mar=c(1,1,0.5,0),new=F,cex=3,lwd=0.5,cex.axis=1,bty="o",mgp=c(3,1.5,0))
# par(mfrow=c(3,2),oma=c(1,1,0.5,1),mar=c(1,1,0.5,1),new=F,cex=1,lwd=0.5,cex.axis=1,bty="o",mgp=c(3,1.5,0))
names(all.ovc)
# Schotia latifolia
# Individual 1
plot(0,col="grey85",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Schotia" & all.ovc$Individual == 1],all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 1],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# 
#Individual 2
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 2],all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 2],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 3
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 4],all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 4],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 4
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 5],all.ovc$PercentEmbolism[all.ovc$Genus  == "Schotia"& all.ovc$Individual == 5],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")

# Schotia
i = 1
# leaves
a <- summary.data[summary.data$tip.label  == "Schotia_latifolia",5]
b <- summary.data[summary.data$tip.label  == "Schotia_latifolia",3]
x <- seq(from=-12,to=0,by=0.05)
y <- 100-100/(1+exp(a*(x-b)))
y.1 <- y
# bootstrapping confidence intervals
pred2 <- function(x,a,b){f <-  100 - 100/(1+exp(a*(x-b)))}
pred.sig <- pred2(x,a,b)
L <- 1000
n.new <- length(x)
df.mc <- data.frame(x=x, f=y)
Y <- matrix(nrow=L,ncol=n.new)
sigmaa.est <- summary.data[summary.data$tip.label  == "Schotia_latifolia",6]*sqrt(summary.data[summary.data$tip.label  == "Schotia_latifolia",7])
sigmab.est <- summary.data[summary.data$tip.label  == "Schotia_latifolia",4]*sqrt(summary.data[summary.data$tip.label  == "Schotia_latifolia",7])
for (l in (1:L)) {
  a.new <- a + rnorm(1,0,sigmaa.est)
  b.new <- b + rnorm(1,0,sigmab.est)
  f.l <- pred2(x,a.new,b.new)
  Y[l,] <- f.l
}
alpha <- 0.05
df.mc[c("lwr.conf","upr.conf")] <- t(apply(Y,MARGIN=2,function(x) quantile(x,c(alpha/2,1-alpha/2))))

x.m.1 <- x
b.p.1 <- b
y.p.1 <- 50
yyy <- y.1[c(max(which(y.1 >99.95)):length(y.1))]
xxx <- x.m.1[c(max(which(y.1 >99.95)):length(y.1))]

# lines(y.2[c(max(which(y.2 >99.95)):length(y.2))]~x.m.2[c(max(which(y.2 >99.95)):length(y.2))],lwd=2,col="lightgrey",lty=1)
lines(yyy~xxx,lwd=1,col="black",lty=1)
points(y.p.1~b.p.1,pch=21,col="black",cex=3,bg=cols[i],ylab="",xlab="")
axis(side=2,las=1,labels=T, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)
axis(side=1,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)

rect(-4.9,-0.075,-5.5,-0.02,col = "grey100",border=NA)
text(-2.5,95,"Schotia latifolia",font=3,cex=1.65)

text(-9.5,90,"a",cex=2,font=1,las=1)
abline(h=0,lty=1,col="black",lwd=0.5)
# tlp data
points(x=tlp.sp[5],y=-3,pch=22,bg="lightgrey",col="black",cex=2)
arrows(tlp.all.n[5],-2.5,tlp.all.p[5],-2.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data 2022
points(x=wp.sp[6],y=-6.5,pch=22,bg="lightblue",col="black",cex=2)
arrows(wp.all.n[6],-6.5,wp.all.p[6],-6.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data 2020
points(x=wp.d.sp[4],y=-10.5,pch=22,bg="orange",col="black",cex=2)
arrows(wp.d.all.n[4],-10.5,wp.d.all.p[4],-10.5,length = 0,code = 3,lwd = 0.5, col="black")
legend(-9.8,40,c("TLP", expression(Psi[min_2022]),expression(Psi[min_2020])),pch=c(22,22,22), col=c("black"),pt.bg = c("lightgrey","lightblue","orange"),cex=1.5) 

#paste(expression(Psi), [min], " 2022")),)
# Polygala myrtifolia
# Individual 1
plot(0,col="grey85",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Polygala" & all.ovc$Individual == 1],all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 1],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 2
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 2],all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 2],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 3
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 3],all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 3],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 4
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 4],all.ovc$PercentEmbolism[all.ovc$Genus  == "Polygala"& all.ovc$Individual == 4],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# Polygala
i <- 2
# leaves
a <- summary.data[summary.data$tip.label  == "Polygala_myrtifolia",5]
b <- summary.data[summary.data$tip.label  == "Polygala_myrtifolia",3]
x <- seq(from=-12,to=0,by=0.05)
y <- 100-100/(1+exp(a*(x-b)))
y.1 <- y
# bootstrapping confidence intervals
pred2 <- function(x,a,b){f <-  100 - 100/(1+exp(a*(x-b)))}
pred.sig <- pred2(x,a,b)
L <- 1000
n.new <- length(x)
df.mc <- data.frame(x=x, f=y)
Y <- matrix(nrow=L,ncol=n.new)
sigmaa.est <- summary.data[summary.data$tip.label  == "Polygala_myrtifolia",6]*sqrt(summary.data[summary.data$tip.label %in% "Polygala_myrtifolia",7])
sigmab.est <- summary.data[summary.data$tip.label  == "Polygala_myrtifolia",4]*sqrt(summary.data[summary.data$tip.label %in% "Polygala_myrtifolia",7])
for (l in (1:L)) {
  a.new <- a + rnorm(1,0,sigmaa.est)
  b.new <- b + rnorm(1,0,sigmab.est)
  f.l <- pred2(x,a,b.new)
  Y[l,] <- f.l
}
alpha <- 0.05
df.mc[c("lwr.conf","upr.conf")] <- t(apply(Y,MARGIN=2,function(x) quantile(x,c(alpha/2,1-alpha/2))))

x.m.1 <- x
b.p.1 <- b
y.p.1 <- 50
yyy <- y.1[c(max(which(y.1 >99.5)):length(y.1))]
xxx <- x.m.1[c(max(which(y.1 >99.5)):length(y.1))]

lines(yyy~xxx,lwd=1,col="black",lty=1)
points(y.p.1~b.p.1,pch=21,col="black",cex=3,bg=cols[i],ylab="",xlab="")
axis(side=2,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)
axis(side=1,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)

rect(-4.9,-0.075,-5.5,-0.02,col = "grey100",border=NA)
text(-2.5,95,"Polygala myrtifolia",font=3,cex=1.65)
abline(h=0,lty=1,col="black",lwd=0.5)
text(-9.5,90,"b",cex=2,font=1,las=1)
# tlp data
points(x=tlp.sp[4],y=-3,pch=22,bg="lightgrey",col="black",cex=2)
arrows(tlp.all.n[4],-2.5,tlp.all.p[4],-2.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.sp[4],y=-6.5,pch=22,bg="lightblue",col="black",cex=2)
arrows(wp.all.n[4],-6.5,wp.all.p[4],-6.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.d.sp[3],y=-10.5,pch=22,bg="orange",col="black",cex=2)
arrows(wp.d.all.n[3],-10.5,wp.d.all.p[3],-10.5,length = 0,code = 3,lwd = 0.5, col="black")

# Pappea capensis
# Individual 1
plot(0,col="grey85",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
par(new=T)
plot(all.ovc$WP[all.ovc$Genus   == "Pappea" & all.ovc$Individual == 2],all.ovc$PercentEmbolism[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 2],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 2
par(new=T)
plot(all.ovc$WP[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 3],all.ovc$PercentEmbolism[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 3],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 3
par(new=T)
plot(all.ovc$WP[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 4],all.ovc$PercentEmbolism[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 4],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 4
par(new=T)
plot(all.ovc$WP[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 5],all.ovc$PercentEmbolism[all.ovc$Genus   == "Pappea"& all.ovc$Individual == 5],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# Pappea
i <- 3
head(summary.data)
# leaves
a <- summary.data[summary.data$tip.label   == "Pappea_capensis",5]
b <- summary.data[summary.data$tip.label   == "Pappea_capensis",3]
x <- seq(from=-12,to=0,by=0.05)
y <- 100-100/(1+exp(a*(x-b)))
y.1 <- y
# bootstrapping confidence intervals
pred2 <- function(x,a,b){f <-  100 - 100/(1+exp(a*(x-b)))}
pred.sig <- pred2(x,a,b)
L <- 1000
n.new <- length(x)
df.mc <- data.frame(x=x, f=y)
Y <- matrix(nrow=L,ncol=n.new)
sigmaa.est <- summary.data[summary.data$tip.label   == "Pappea_capensis",6]*sqrt(summary.data[summary.data$tip.label   == "Pappea_capensis",7])
sigmab.est <- summary.data[summary.data$tip.label   == "Pappea_capensis",4]*sqrt(summary.data[summary.data$tip.label   == "Pappea_capensis",7])
for (l in (1:L)) {
  a.new <- a + rnorm(1,0,sigmaa.est)
  b.new <- b + rnorm(1,0,sigmab.est)
  f.l <- pred2(x,a.new,b.new)
  Y[l,] <- f.l
}
alpha <- 0.05
df.mc[c("lwr.conf","upr.conf")] <- t(apply(Y,MARGIN=2,function(x) quantile(x,c(alpha/2,1-alpha/2))))

x.m.1 <- x
b.p.1 <- b
y.p.1 <- 50
yyy <- y.1[c(max(which(y.1 >99.95)):length(y.1))]
xxx <- x.m.1[c(max(which(y.1 >99.95)):length(y.1))]
lines(yyy~xxx,lwd=1,col="black",lty=1)
points(y.p.1~b.p.1,pch=21,col="black",cex=3,bg=cols[i],ylab="",xlab="")
axis(side=2,las=1,labels=T, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)
axis(side=1,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)

rect(-4.9,-0.075,-5.5,-0.02,col = "grey100",border=NA)
text(-2.5,95,"Pappea capensis",font=3,cex=1.65)
abline(h=0,lty=1,col="black",lwd=0.5)
text(-9.5,90,"c",cex=2,font=1,las=1)
# tlp data
points(x=tlp.sp[3],y=-3,pch=22,bg="lightgrey",col="black",cex=2)
arrows(tlp.all.n[3],-2.5,tlp.all.p[3],-2.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.sp[3],y=-6.5,pch=22,bg="lightblue",col="black",cex=2)
arrows(wp.all.n[3],-6.5,wp.all.p[3],-6.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.d.sp[2],y=-10.5,pch=22,bg="orange",col="black",cex=2)
arrows(wp.d.all.n[2],-10.5,wp.d.all.p[2],-10.5,length = 0,code = 3,lwd = 0.5, col="black")

# Boscia oleoides
par(new=F)
i <- 3
plot(0,col="grey85",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# Individual 1
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Boscia" & all.ovc$Individual == 1],all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 1],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 2
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 2],all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 2],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 3
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 3],all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 3],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 4
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 4],all.ovc$PercentEmbolism[all.ovc$Genus  == "Boscia"& all.ovc$Individual == 4],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# Boscia
# # mean curve
a <- summary.data[summary.data$tip.label  == "Boscia_oleoides",5]
b <- summary.data[summary.data$tip.label  == "Boscia_oleoides",3]
x <- seq(from=-12,to=0,by=0.05)
y <- 100-100/(1+exp(a*(x-b)))
y.1 <- y
# bootstrapping confidence intervals
pred2 <- function(x,a,b){f <-  100 - 100/(1+exp(a*(x-b)))}
pred.sig <- pred2(x,a,b)
L <- 1000
n.new <- length(x)
df.mc <- data.frame(x=x, f=y)
Y <- matrix(nrow=L,ncol=n.new)
sigmaa.est <- summary.data[summary.data$tip.label  == "Boscia_oleoides",6]*sqrt(summary.data[summary.data$tip.label  == "Boscia_oleoides",7])
sigmab.est <- summary.data[summary.data$tip.label  == "Boscia_oleoides",4]*sqrt(summary.data[summary.data$tip.label  == "Boscia_oleoides",7])
for (l in (1:L)) {
  a.new <- a + rnorm(1,0,sigmaa.est)
  b.new <- b + rnorm(1,0,sigmab.est)
  f.l <- pred2(x,a.new,b.new)
  Y[l,] <- f.l
}
alpha <- 0.05
df.mc[c("lwr.conf","upr.conf")] <- t(apply(Y,MARGIN=2,function(x) quantile(x,c(alpha/2,1-alpha/2))))

x.m.1 <- x
b.p.1 <- b
y.p.1 <- 50
yyy <- y.1[c(max(which(y.1 >99.95)):length(y.1))]
xxx <- x.m.1[c(max(which(y.1 >99.95)):length(y.1))]
lines(yyy~xxx,lwd=1,col="black",lty=1)
points(y.p.1~b.p.1,pch=21,col="black",cex=3,bg=cols[i])
axis(side=2,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)
axis(side=1,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)

rect(-4.9,-0.065,-6,-0.02,col = "grey100",border=NA)
text(-2.5,95,"Boscia oleoides",font=3,cex=1.65)
abline(h=0,lty=1,col="black",lwd=0.5)
text(-9.5,90,"d",cex=2,font=1,las=1)
# tlp data
points(x=tlp.sp[1],y=-3,pch=22,bg="lightgrey",col="black",cex=2)
arrows(tlp.all.n[1],-2.5,tlp.all.p[1],-2.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.sp[1],y=-6.5,pch=22,bg="lightblue",col="black",cex=2)
arrows(wp.all.n[1],-6.5,wp.all.p[1],-6.5,length = 0,code = 3,lwd = 0.5, col="black")

# Euclea undulata
# Individual 1
plot(0,col="grey85",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Euclea" & all.ovc$Individual == 1],all.ovc$PercentEmbolism[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 1],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 2
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 2],all.ovc$PercentEmbolism[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 2],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 3
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 3],all.ovc$PercentEmbolism[all.ovc$Genus  == "Euclea"& all.ovc$Individual == 3],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# Euclea
i <- 5
# leaves
a <- summary.data[summary.data$tip.label  == "Euclea_undulata",5]
b <- summary.data[summary.data$tip.label  == "Euclea_undulata",3]
x <- seq(from=-12,to=0,by=0.05)
y <- 100-100/(1+exp(a*(x-b)))
y.1 <- y
# bootstrapping confidence intervals
pred2 <- function(x,a,b){f <-  100 - 100/(1+exp(a*(x-b)))}
pred.sig <- pred2(x,a,b)
L <- 1000
n.new <- length(x)
df.mc <- data.frame(x=x, f=y)
Y <- matrix(nrow=L,ncol=n.new)
sigmaa.est <- summary.data[summary.data$tip.label  == "Euclea_undulata",6]*sqrt(summary.data[summary.data$tip.label  == "Euclea_undulata",7])
sigmab.est <- summary.data[summary.data$tip.label  == "Euclea_undulata",4]*sqrt(summary.data[summary.data$tip.label  == "Euclea_undulata",7])
for (l in (1:L)) {
  a.new <- a + rnorm(1,0,sigmaa.est)
  b.new <- b + rnorm(1,0,sigmab.est)
  f.l <- pred2(x,a.new,b.new)
  Y[l,] <- f.l
}
alpha <- 0.05
df.mc[c("lwr.conf","upr.conf")] <- t(apply(Y,MARGIN=2,function(x) quantile(x,c(alpha/2,1-alpha/2))))

x.m.1 <- x
b.p.1 <- b
y.p.1 <- 50
yyy <- y.1[c(max(which(y.1 >99.95)):length(y.1))]
xxx <- x.m.1[c(max(which(y.1 >99.95)):length(y.1))]

# lines(y.2[c(max(which(y.2 >99.95)):length(y.2))]~x.m.2[c(max(which(y.2 >99.95)):length(y.2))],lwd=2,col="lightgrey",lty=1)
lines(yyy~xxx,lwd=1,col="black",lty=1)
points(y.p.1~b.p.1,pch=21,col="black",cex=3,bg=cols[i],ylab="",xlab="")
axis(side=2,las=1,labels=T, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)
axis(side=1,las=1,labels=T, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)

rect(-4.9,-0.075,-5.5,-0.02,col = "grey100",border=NA)
text(-2.5,95,"Euclea undulata",font=3,cex=1.65)
abline(h=0,lty=1,col="black",lwd=0.5)
text(-9.5,90,"e",cex=2,font=1,las=1)
# tlp data
points(x=tlp.sp[2],y=-3,pch=22,bg="lightgrey",col="black",cex=2)
arrows(tlp.all.n[2],-2.5,tlp.all.p[2],-2.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.sp[2],y=-6.5,pch=22,bg="lightblue",col="black",cex=2)
arrows(wp.all.n[2],-6.5,wp.all.p[2],-6.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.d.sp[1],y=-10.5,pch=22,bg="orange",col="black",cex=2)
arrows(wp.d.all.n[1],-10.5,wp.d.all.p[1],-10.5,length = 0,code = 3,lwd = 0.5, col="black")

# Searsia longispina
par(new=F)
i <- 6
plot(0,col="grey85",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# Individual 1
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Searsia" & all.ovc$Individual == 1],all.ovc$PercentEmbolism[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 1],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 2
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 2],all.ovc$PercentEmbolism[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 2],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
#Individual 3
par(new=T)
plot(all.ovc$WP[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 3],all.ovc$PercentEmbolism[all.ovc$Genus  == "Searsia"& all.ovc$Individual == 3],col="grey50",yaxt="n",xlim=c(-10,0),ylim=c(-12,100),ylab="",xlab="",xaxt="n",type="l")
# mean curve
a <- summary.data[summary.data$tip.label  == "Searsia_longispina",5]
b <- summary.data[summary.data$tip.label  == "Searsia_longispina",3]
x <- seq(from=-12,to=0,by=0.05)
y <- 100-100/(1+exp(a*(x-b)))
y.1 <- y
# bootstrapping confidence intervals
pred2 <- function(x,a,b){f <-  100 - 100/(1+exp(a*(x-b)))}
pred.sig <- pred2(x,a,b)
L <- 1000
n.new <- length(x)
df.mc <- data.frame(x=x, f=y)
Y <- matrix(nrow=L,ncol=n.new)
sigmaa.est <- summary.data[summary.data$tip.label  == "Searsia_longispina",6]*sqrt(summary.data[summary.data$tip.label  == "Searsia_longispina",7])
sigmab.est <- summary.data[summary.data$tip.label  == "Searsia_longispina",4]*sqrt(summary.data[summary.data$tip.label  == "Searsia_longispina",7])
for (l in (1:L)) {
  a.new <- a + rnorm(1,0,sigmaa.est)
  b.new <- b + rnorm(1,0,sigmab.est)
  f.l <- pred2(x,a.new,b.new)
  Y[l,] <- f.l
}
alpha <- 0.05
df.mc[c("lwr.conf","upr.conf")] <- t(apply(Y,MARGIN=2,function(x) quantile(x,c(alpha/2,1-alpha/2))))

x.m.1 <- x
b.p.1 <- b
y.p.1 <- 50
yyy <- y.1[c(max(which(y.1 >99.95)):length(y.1))]
xxx <- x.m.1[c(max(which(y.1 >99.95)):length(y.1))]
lines(yyy~xxx,lwd=1,col="black",lty=1)
points(y.p.1~b.p.1,pch=21,col="black",cex=3,bg=cols[i])
axis(side=2,las=1,labels=F, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)
axis(side=1,las=1,labels=T, tick=T,lwd=0.5,tck=-0.025,cex.axis=1.5)

rect(-4.9,-0.065,-6,-0.02,col = "grey100",border=NA)
text(-2.5,95,"Searsia longispina",font=3,cex=1.65)
abline(h=0,lty=1,col="black",lwd=0.5)
text(-9.5,90,"f",cex=2,font=1,las=1)
# tlp data
points(x=tlp.sp[6],y=-3,pch=22,bg="lightgrey",col="black",cex=2)
arrows(tlp.all.n[6],-2.5,tlp.all.p[6],-2.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.sp[7],y=-6.5,pch=22,bg="lightblue",col="black",cex=2)
arrows(wp.all.n[7],-6.5,wp.all.p[7],-6.5,length = 0,code = 3,lwd = 0.5, col="black")
# water potential data
points(x=wp.d.sp[5],y=-10.5,pch=22,bg="orange",col="black",cex=2)
arrows(wp.d.all.n[5],-10.5,wp.d.all.p[5],-10.5,length = 0,code = 3,lwd = 0.5, col="black")

mtext(side=2,"Observed embolism (%)",outer=T,line=3,at=0.5,cex=7)
mtext(side=1,"Xylem water potential (MPa)",outer=T,line=3.75,at=0.5,cex=7)
dev.off()